knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
#devtools::install_github('IQSS/Zelig')
library(readstata13)
library(Zelig)
library(stargazer)
library(latex2exp)
library(dplyr)
#ref for relogit: code: http://docs.zeligproject.org/articles/zelig_relogit.html
#explanation: https://www.r-bloggers.com/2018/06/bias-adjustment-for-rare-events-logistic-regression-in-r/

directory = c("~/Dropbox (Harvard University)/Coursework/Harvard/2022 Fall/Gov 2001/Climate conflict")
rep_df <- read.dta13(paste0(directory[1],"/replication materials/Nel and Righarts ND and VCC final 20 June 2007.dta"))
rep_df <- rep_df %>%                           
  group_by(state) %>%
  mutate(imr_t_1 = lag(imr),
         imr_sq_t_1 = lag(imr_sq),
         vccall_tp1 = lead(vccall),
         vccmajor_tp1 = lead(vccmajor)) 

## Table 1
#Model 1
t1_1 <- zelig(vccall ~ allND_pc + imr_t_1 + imr_sq_t_1 + mixedx + gdpgrox + brevity, 
              model = "relogit", data = rep_df)

#Model 2
t1_2 <- zelig(vccall ~ allND_pc + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
              model = "relogit", data = rep_df)

#Model 3
t1_3 <- zelig(vccall ~ allND_pc + gdp_pc + mixed + gdpgro + brevity, 
              model = "relogit", data = rep_df)

#Model 4
t1_4 <- zelig(vccall ~ allND_pc + mixed + gdpgro + youth_bulge + brevity, 
              model = "relogit", data = rep_df)

#Model 5
t1_5 <- zelig(vccall_tp1 ~ allND_pc + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
              model = "relogit", data = rep_df)


stargazer(from_zelig_model(t1_1), from_zelig_model(t1_2), from_zelig_model(t1_3), from_zelig_model(t1_4), from_zelig_model(t1_5), title="Natural Disasters and the Risk of Violent Civil Conflict, 1950-2000", column.labels=c("\\specialcell{All violent\\\\civil conflict}","\\specialcell{All violent\\\\civil conflict}","\\specialcell{All violent\\\\civil conflict}","\\specialcell{All violent\\\\civil conflict}","\\specialcell{All violent\\\\civil conflict $t+1$}"),
          covariate.labels=c("No. natural disasters per capita","Infant mortality rate $t-1$",
                             "Infant mortality rate squared $t-1$", "Mixed regime (m)","GDP growth (m)","GDP per capita","Mixed regime","GDP growth","Youth bulge", "Brevity of peace"), dep.var.labels.include = FALSE, no.space=TRUE,
          omit.stat=c("LL","ser","f","aic"), align=TRUE, out = paste0(directory[1],"/output/tables/rep_table1.tex"), header=FALSE )


## Table2
#Model 6
t2_6 <- zelig(vccall ~ geolND + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
              model = "relogit", data = rep_df)

#Model 7
t2_7 <- zelig(vccall ~ geolND + mixed + gdpgro + youth_bulge + brevity + size_ln, 
              model = "relogit", data = rep_df)

#Model 8
t2_8 <- zelig(vccminor ~ geolND + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
              model = "relogit", data = rep_df)

#Model 9
t2_9 <- zelig(vccmajor ~ geolND + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
              model = "relogit", data = rep_df)

#Model 10 
t2_10 <- zelig(vccmajor_tp1 ~ geolND + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
               model = "relogit", data = rep_df)

stargazer(from_zelig_model(t2_6), from_zelig_model(t2_7), from_zelig_model(t2_8), from_zelig_model(t2_9), from_zelig_model(t2_10), title="Geological Disasters and the Risk of Violent Civil Conflict, 1950-2000", column.labels=c("\\specialcell{All violent\\\\civil conflict}","\\specialcell{All violent\\\\civil conflict}","\\specialcell{Minor violent\\\\civil conflict}","\\specialcell{Major violent\\\\civil conflict}","\\specialcell{Major violent\\\\civil conflict $t+1$}"),
          covariate.labels=c("Incidence of geological disasters","Infant mortality rate $t-1$",
                             "Infant mortality rate squared $t-1$","Mixed regime","GDP growth","Youth bulge", "Brevity of peace", "Size"), 
          dep.var.labels.include = FALSE, no.space=TRUE,
          omit.stat=c("LL","ser","f","aic"), align=TRUE, out = paste0(directory[1],"/output/tables/rep_table2.tex"), header=FALSE)

## Table 3
#Model 11
t3_11 <- zelig(vccall ~ allclim_pc + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
               model = "relogit", data = rep_df)

#Model 12
t3_12 <- zelig(vccminor ~ allclim_pc + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
               model = "relogit", data = rep_df)

#Model 13
t3_13 <- zelig(vccmajor ~ allclim_pc + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
               model = "relogit", data = rep_df)

#Model 14
t3_14 <- zelig(vccmajor_tp1 ~ allclim_pc + imr_t_1 + imr_sq_t_1 + mixed + gdpgro + brevity, 
               model = "relogit", data = rep_df)

stargazer(from_zelig_model(t3_11), from_zelig_model(t3_12), from_zelig_model(t3_13), from_zelig_model(t3_14), title="Climate-related Disasters and the Risk of Violent Civil Conflict, 1950-2000", column.labels=c("\\specialcell{All violent\\\\civil conflict}","\\specialcell{Minor violent\\\\civil conflict}","\\specialcell{Major violent\\\\civil conflict}","\\specialcell{Major violent\\\\civil conflict $t+1$}"),
          covariate.labels=c("No. climate-related natural disasters per capita","Infant mortality rate $t-1$",
                             "Infant mortality rate squared $t-1$","Mixed regime","GDP growth", "Brevity of peace"), 
          dep.var.labels.include = FALSE, no.space=TRUE,
          omit.stat=c("LL","ser","f","aic"), align=TRUE, out = paste0(directory[1],"/output/tables/rep_table3.tex"), header=FALSE)
